Late-time evolution of the gravitating Skyrmion 
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We study the dynamics of spherically symmetric solutions in the Einstein-Skyrme model. We 
focus our attention on generic long-time evolution of initial data resulting in the formation of the 
B — 1 soliton, which plays the role of an attractor. We demonstrate that similarly to the case of 
flat space evolution, the relaxation to the regular soliton (which we will call Skyrmion) is universal 
and may be treated as a superposition of two effects - quasinormal oscillations responsible for 
intermediate asymptotics and a power-law tail describing the behavior of the system at very long 
times. We determine the values of parameters describing asymptotics and examine their dependence 
on the value of dimensionless coupling constant of the model. 
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I. INTRODUCTION 
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This paper is concerned with intermediate and final asymptotics of spherically symmetric solutions in Einstein- 
. Skyrme (ES) model and is meant as a natural extension of similar analysis done in ES model is known 0, Q 

■ as the model which admits a surprisingly rich spectrum of asymptotically flat static spherically symmetric solutions, 
both with and without horizon. In particular, the static B = 1 spherically symmetric solution (Skyrmion) is stable 

' [3, 0, Q • Such stable static solutions are very important configurations in any physical model as they are natural 
candidates for endstates of a generic evolution. During this evolution the system relaxes to the final stable configuration 
^ ' radiating energy to infinity. This radiation contains information about the nature of the central object, therefore the 
l^f^ understanding of the mechanism of the emission is of crucial importance. 

'— ' Skyrme model was suggested by T.H.R. Skyrme [6| as a model describing properties of baryons in particle physics. 

^ It is a nonlinear model of a chiral field in which baryons are identified as solitons - topologically stable static field 
^ configurations. Apart from its usefulness in particle physics the model is interesting theoretically because despite 
' of its simplicity we may observe here - both in flat space as well as in self-gravitating case - all features of the 
fS| , relaxation process. So for generic initial data after short period of evolution sensitive to the initial conditions, the 
' system reaches intermediate asymptotics which may be described by the sum of exponentially decaying oscillations 
\l [ (so called quasinormal modes) superimposed on final, static Skyrmion configuration. Its long-time asymptotics is 
' governed by the power-law tail. 

■ The rest of this paper is organized as follows. We start in Section II with the specification of the model. Then in 
0^ . Section III we discuss the linear stability of the Skyrmion. Next Section IV is devoted to the calculation of parameters 

' of quasinormal modes. In Section V we describe the numerical method used in calculation of the time evolution and 
^ . demonstrate the results of calculation of quasinormal modes obtained directly from the numerical evolution as well 
as from the method described in Section IV. In this Section we discuss also the results for power-law tails. 



II. SETTINGS 



We consider the Einstein-Skyrme model, so the matter in our model is a chiral field - an SU{2) - valued scalar 
function U{x) with dynamics given by the Lagrangian Q: 

L = trriV^V'^U-') + ^Tr[{VaU)U-\ (V,C/)f/-T - j^R, (1) 

where Vq is the covariant derivative with respect to the spacetime metric, G - gravitational constant and R - scalar 
of curvature. By choosing appropriate units of mass and length we set the values of coupling constants / and e to 
one. 

We restrict ourselves to the case of spherical symmetry and use the polar time slicing and the areal radial coordinate 
so we may parameterize the metric as follows 

ds'^ = -e-2''('''*)iV(r, t)dt^ + N-\r, t)dr^ + r^dQ^, (2) 

where dfl^ is the standard metric on the unit 2-sphere. For the chiral field we apply the usual hedgehog ansatz 
U = exp(i~a ■ fF{r, t)), where It is the vector of Pauli matrices and f - unit radial vector. 
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We denote ^ and ^ by overdots and primes respectively and introduce an auxiliary variable P = ue^N ^F, where 
M = + 2 sin^ F and as a result obtain the following set of ES equations: 

F^e~^N-, (3) 
u 



P = {e-'NuFj + sin(2^^)e-^ Ar(^ - F'^) . 1 , (4) 



N = -—e-^N^PF', (5) 



r 



S'--^(^ + F-], (6) 



r \u 

N' = l^^^(2 sin^ F + ^ + .iV(^ + F'^)) . (7) 

Here a = AttGJ^ is dimensionless coupling constant. 

We are interested in regular asymptotically flat solution of ES equations. To ensure the regularity in the center, 
N{r) = 1 + O(r^), we impose the following boundary condition F{r,t) ^ r for r — > 0. Asymptotic flatness of 
initial data N{r,0) = 1 + 0{l/r) requires that F{r,0) = Bit + 0(l/r^) at infinity, where the integer B, which we 
will call baryon number, is equal to the topological degree of the chiral fleld. As long as no horizon forms the baryon 
number of initial configuration is preserved during the evolution. Therefore the asymptotic flatness condition breaks 
the initial value problems into infinitely many disjoint topological sectors labeled by the baryon number B. Below we 
consider mainly B = 1 sector. In is a well established fact |3j that in this sector for a < acrit — 0.040378 the Eqs.(l3][7|) 
have two regular static solutions. One of them is linearly stable and is identical with flat Skyrmion for a = 0, so we 
regard it just as a gravitationally distorted Skyrmion. The second one, which we denote by X" is unstable. This two 
solutions coalesce at a — acrit and disappear for a > acrit- 



III. LINEAR STABILITY OF GRAVITATING SKYRMION 



We are interested in relaxation processes of solutions to the static configuration, therefore we start with the linear 
stability analysis of Skyrmion. As we restrict ourselves only to the spherically symmetric perturbations, we use the 
following standard decomposition of solutions: 

F{r,t) ^ S{r) + h{r,t), 7V(r, t) = A^o W + i), 6{r,t) ^ So{r) + Si{r,t), (8) 

where (S'(r), Ao(r), (5o(f )) denote static Skyrmion configuration. We plug this decomposition into Eqs.([3][7]) and 
linearize. As it was demonstrated in the analysis of such system simplifies considerably because in spherical 
symmetry the metric perturbations which enter the pulsation equation for /i (r, t) are completely determined by the 
matter perturbations. Exploiting this and introducing auxiliary field v{r,t) defined by the formula 

^(M) = ^=^^M=, (9) 
Vr2+ 2 sinks' 

we obtain the linear pulsation equation for u(r, t): 

e^°N-'^v- {e-^"Nov'y+ e'^^Uv = 0. (10) 

Let us note that the potential U — Uq + aUi + is determined entirely by the functions describing the static 

Skyrmion and is given by the formulae: 
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Ui^8w{l + w^)S'cosS- '^l^l\^2^ - 2(1+ 2w^)S'^ (12) 
U2^2{1+ w^) (2+ w^)S''^ sin^ S. (13) 

Here w = sin(^). 

Let us examine the behavior of the finite part of potential V{r) = U{r) — 2/r^ for large r. We may easily find, 
that this behavior in the fiat space and in the case of the self-gravitating model are different. In both cases we have 
S{r) ^ TT — 1/r^ what means that w ^ In case of the flat model we have additionally a = 0, NQ{r) — \\ this 

means that potential V{r) ~ l/r^. For self-gravitating case A'o ~ 1 — h/r what leads to the relation V{r) ^ 1/r^ for 
r — > 0. This difference will be important in the discussion of power-law tails. 

Using the radial coordinate p defined by 

^ = e^«W-^«(°)iVo-i, p(0) = 0, (14) 
dr 

and assuming the time-dependence of the form i;(r, t) = e~*'^*'0(p)7 we transfer Eq. (llO|) into a radial p-wave Schrodinger 
equation ^] 

(-^ + ^+ np)) = feV^-^^^V-. (15) 



F(p)=e-2^"M+2^°(")f/(r(p))- ^. (16) 



Here the potential 



is bounded for all p. The investigation of linear stability is then reduced to study the eigenvalue problem (jlSp in 
the space of square integrable functions. This was done in where it was found, that the spectrum around the 
Skyrmion is continuous and positive, > 0. This means that the gravitating Skyrmion is linearly stable. If we make 
similar calculation for X" we find out that Eq. p^ has exactly one bound state with A:^ < 0, indicating an instability 
- exponentially growing mode with the exponent = ikb > 0. In addition, for a (Xcrit, the instability exponent 



le exp 

7 — > 0. Numerical calculation done in Q confirm the above results also in nonlinear analysis. 



IV. CALCULATION OF QUASINORMAL MODES 



Linear stability of the Skyrmion is an important property of ES system because it tells us that this solution is an 
attractor describing the final configurations in _B = 1 sector of initial data problem. However it does not provide a 
direct information about the way this configuration is reached in the evolution. To obtain the intermediate and long 
time asymptotics of solutions of Eqs.(l3][7l) one usually discuss the Eq. p^ describing the perturbation of Skyrmion. 
A problem of this type - especially perturbations of black holes and compact objects - has a long history starting 
with the papers of Regge and Wheeler 9], Zerilli [10] and Price [Tl[, for a review see [1] and [l3- As we will see the 
perturbations of Skyrmion at intermediate times are described by quasinormal modes, whereas late time asymptotics 
is controlled by a power-law tail. To simplify the notation we will neglect the difference between the radial coordinates 
r and p, so the Eq. pS]) takes the form: 



*" + (4 + V(r))*^fc2^. (^7) 



In the case of ES system by quasinormal mode we will understand the solution of Eq. (|17p satisfying the outgoing 
wave conditions for r oo: 



*(r) oc e^'^'', k = u>-ij, 7 > 0, (18) 

and regular at the center. 

We expect, that the mode with the least damping 7 will dominate the intermediate asymptotics in the B = 1 
sector. To find a proper solution we have to apply the above boundary condition to the solution of Eq. (|17p : this will 
quantize the eigenvalue k. 
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However, it is not easy to implement the boundary condition of this type (especially for long-range potentials). To 
see that let us remark, that for r > R, where R - the range of potential, any regular solutions 5" of Eq. pT|) may be 

decomposed in terms of Riccati-Hankel functions h\^\kr) and h\ \kr) as follows: 

* = a+{k) }i-^\kr) + a_(fc) }ii\kr). (19) 

The Riccati-Hankel functions have the following asymptotic behavior: hf{kr) ^ e^*'^'', hence the outgoing wave 
boundary conditions correspond to the zeros of a-{k) coefficient. However, this is numerically difficult to achieve as 
the "unwanted" ingoing component of (|19p decreases exponentially with r. Therefore, we try to remove a numerically 
small ingoing component on the background of large outgoing component. To achieve this we should have there a very 
good resolution, at least of order of 0(e~^^^). For this reason a usual procedure of shooting from r = and varying 
k until a-(k) vanishes (naive shooting) is not feasible, especially in case of ^R ^ 1. We have found that much better 
way of proceeding is to use a shooting-to-a-fitting-point technique which is much more accurate and additionally may 
be generalized to the case of long-range potentials. In this technique we shoot from two sides - r = (requiring that 
the solution is regular for r — > 0) and some large r2, where the decomposition (jl9p works and try to match both 
solutions at some intermediate rf. 

Technically [H, we use the amplitude-phase representation for the solution of Eq. p?)) : "if — A exp{i^). As a result, 
this equation for complex function \1/ is replaced by the following system of two equations for real functions A and $: 

A"-A(j,'^~(^ + V + -f^-LjAA = 0, (20) 



Regularity at the center requires: 



A(j)" + 2A'0' - 2(jjjA = 0. (21) 



A{r) - and <j){r) - —r^ for r ^ 0. (22) 
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We solve this equations numerically from r = up to some relatively small intermediate rf. For r > rf we replace 
the Schrodinger equation (fT7|) by its Riccati form 

g' + - 2/r^ -V + k^ = 0, (23) 

where g is defined as a logarithmic derivative of ^': 3 = ^ = ^ + We solve this equation backwards in r from 

some large r2 to the intermediate rf defined above, starting with the initial value g{r2) — ^ • We assume, 

that r2 is large enough, so the Riccati-Hankel function h'l^\kr) — (— i + -^)e^^^ which is exact outgoing solution of 
the free {V — 0) Riccati equation approximates well the solution of the full Riccati equation. Now the quantization 
condition for the quasinormal modes is the matching of logarithmic derivatives at intermediate matching point r / . 



V. NUMERICS AND RESULTS 



We have performed extensive numerical studies of long time asymptotics of solutions in ES model. To this end we 
have solved initial value problem ([3][7]) for various baryon numbers, coupling constant and initial data. To do this we 
apply the standard method of lines to the pair of dynamic equations ([3]) and This means that we discretize space 
and replace spatial derivatives by proper algebraic difference approximations. In that way we change the original 
set of dynamical PDEs into a system of ODEs, which may be solved by standard methods. Usually we use 5-point, 
fourth-order accurate spatial discretization and solve the resulting ODE system by means of fourth-order Runge-Kutta 
method. In addition, on each time layer we update the metric functions N and 5 by solving the slicing condition ^ 
and hamiltonian constraint ([7]). Here we also use fourth-order Runge-Kutta method with spline interpolation for the 
values of grid functions at the positions out of the grid which are also required by the integration procedure. As a 
result we have finite difference method which is fourth-order accurate both in space and time. We apply this method 
to the calculation of time evolution which starts with different initial conditions. Typical examples of initial data are: 

F(i = 0,r) = tanh(a;/s), P(t = 0,r) = (24) 

for i? = 1 sector and 

F{t^Q,r)=Ar^cyi^[-{{r-rQ)/s)% P{t = Q,r)^Q (25) 
for B = sector of initial value problem. 
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To ensure the regularity at the origin we impose the boundary conditions: F{t, r = 0) = 0, P{t, r = 0) = 0. On the 
outer boundary we apply outgoing wave boundary condition. This is not enough when we calculate the properties of 
tails; in order to avoid the contamination of results by the part of solution reflected from the outer boundary, we use 
the size of the grid big enough, so the calculation stops before the reflected signal reaches the observation point. In 
calculations of this kind the quadruple precision is crucial, otherwise the accumulation of round-off errors spoils the 
results at late times. 

The typical results of long time evolution in B — 1 sector is demonstrated in Fig[T]a. In this figure a value of the 
field P measured at a fixed value of r is plotted as a function of t. In this figure we may observed all features of 
relaxation process towards a stable Skyrmion configuration: the beginning of evolution is strongly dependent on the 
shape of initial data. After some time an intermediate asymptotics sets in (5 < t < 100) where the relaxation manifest 
itself as a damped oscillations. Finally, at large times {t > 100) we observe the power-law tail. The same features 
may be observed on a complementary picture where a late-time snapshot of solution (see Fig[TJb), i.e. \P{t = ta,r)\ 
as a function of r is seen. Here we may also observed, that the fragment of solution corresponding to the quasi-normal 
modes treated as a function of r grows exponentially, what makes the procedure of finding the mode parameters 
difhcult numerically. Using numerical data we may extract the parameters of least-damped quasi-normal mode (see 
FiglUfor details). 




I ^ ^ , ^ ^ , 1 ie_ig I , ^ ^ ^ ^ ^ ^ , 1 

50 100 150 200 250 300 20 40 60 80 100 120 140 160 180 

t r 



FIG. 1: An example of time evolution in B = 1 sector. Left panel: we plot ln\P{t,ro = 5)| as a function of time and observe 
different phases of relaxation process. Right panel: the snapshot of a solution for ta ~ 160 




FIG. 2: Intermediate asymptotics of relaxation towards Skyrmion. Fitting the exponentially damped oscillation P{t, ro) = 
Ae~'^^\cos{ujt + a)\ to the numerical data we may estimate the parameters of the least damped quasi-normal mode. 



To calculate the parameters of quasi-normal modes we use also the semi-analytic method described in Section IV. 
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To achieve this, we solve - using standard fourth-order Runge-Kutta method - Eqs. ((20)) . ((2T|) . (|23p together with 
time-independent version of system ([2HZI), which gives the potential V. To get the value of (w,7) parameters of the 
least damped quasi-normal mode we proceed as follow: we treat the method described in Section IV as a function T 
which transforms starting values of mode parameters into the values corresponding to the mode in question: 



(26) 



We require that the true solution corresponds to a fixed point of transformation T with the universal values of starting 
mode parameters wq = l;7o = 1- In addition we require that the T procedure is not sensitive to the values of the 
intermediate matching point rf and external shooting point r2. This means, that if we look at the value of the solution 
as a function of and r2 we should observe a plateau. Figl3] shows that it is really the case. We observe, that for 
r2 — 14, the solution of semi-analytic method does not depend on the value of for 1 < < 3. Similarly - if we fix 
rf = 2 than the solution of semi-analytic method does not depend on the value of r2 for 9 < r2 < 18. In addition in 
both cases the value of semi-analytic solution agrees with the solution obtained as a fit to the numerical data. This 
confirms that the semi-analytic method works well and is a good tool for estimation of the parameters of the least 
damped quasi-normal mode. 
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FIG. 3: The sensitivity of results for the parameters (cu,^) obtained via shooting-to-a-fitting-point technique on and r2 
parameters. Flat regions in both plots demonstrate the robustness of the method. 

Values of the QNR parameters as a function of dimensionless coupling constant a are plotted in FigH) Both 
parameters decrease with growing a and go to zero as a approaches the critical value acrit- It means that, as 
a — > acrit the wavelength of perturbation of Skyrmion and life-time of these perturbations are growing to oo. 




FIG. 4: The dependence of the parameters of the fundamental quasinormal mode on the coupling constant. The solid lines 
and points denote data obtained from fit to numerical solutions and by means of shooting technique, respectively. 
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In Figl5]we present the dispersion relation 7 vs. lo for a — > Ucrit- Numerical analysis seems to suggest, that this 
relations takes the form: 7 ^ cj"^. 



-2.5 



-3 - 



-3.5 



-4.5 



fit to data 



log(y) = a*log(B)) + b 
a = 2.9909 
b =-0.5341 



-1.5 -1.4 -1.3 -1.2 -1.1 -1 -0.9 -0.£ 

log(Q)) 



FIG. 5; The dispersion relation at the critical value of coupling constant a — > Ucri 



We have also studied late time asymptotics, where quasi-normal modes are negligible and the relaxation process is 
dominated by power-law tail. Our numerical calculations show that at late times the relaxation to the Skyrmion takes 
the form: F{t^ r) — S{r) t^'^ where the value of 7 depends on the coupling constant a and is equal 5 for a = and 
4 for a > 0. As it was already pointed out in [l[ for the case of flat space, this does not agree with the predictions of 
linear scattering theory (see e.g. [H!, According to this theory, for late times the evolution of the system is well 
described by a hnear wave equation with potential V{r) (see Eq. lfTO]) ). Let us consider more general case, described 
by the equation: 



1/ * = fc^^f. 



(27) 



where V{r) is the finite part of potential. If the fall-off of V{r) at spatial infinity is /3, i.e. V{r) ^ r^'^ for r — *■ 00, 
than for compactly supported initial data linear theory predicts: 7 = 2/ -I- /3. As in our problem we have / = 1 and 
/? = 6 for flat and /3 = 3 for gravitating case (see Section III) , we may expect 8 and 5 for the values of power-law tail 
exponent for flat and gravitating model respectively, which is in disagreement with numerical results. This is another 
example of a situation where the tail has genuinely non-linear character and linear theory fails. To describe the tail 
correctly, one should use a nonlinear perturbatively scheme proposed recently by Bizon et. al. [l|, [l^ [l^, [13, [3 • The 
application of this technique to the case of i? = 1 sector of ES model is tedious, so we have checked, that we get the 
same values of tail exponents in the case of i? = sector of ES model and in non-linear sigma model. Analyzes of 
tails in these models will be published elsewhere (see fl^ for more details). 
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